AD-A080  311 
UNCLASSIFIED 


RENSSELAER  POLYTECHNIC  INST  TROY  NY  DEPT  OF  MATHEMAT— ETC  F/©  12/1 
COLLOCATION  WITH  POLYNOMIAL  AND  TAUT  SPLINES  FOR  SIN6ULARLY  PER— ETC(U> 
NOV  79  J  E  FLAHERTY t  W  MATHON  AFOSR-75-2818 

AFOSR-TR-BO-OOS5  NL 


8JA08031 1 


r 


l.  r:  ■  ,  J 

£OLLOCATION  WITH  POLYNOMIAL  AND  TAUT 
SPLINES  FOR  SINGULARLY  PERTURBED 


BOUNDARY  VALUE  PROBLEMS *r  / 

TT  •  - / 

...  ^ 

\i 0  Joseph  E.^ Flaherty  ‘  ^ 


William 'Mathon 

Department  of  Mathematical  Sciences 
Rensselaer  Polytechnic  Insti 
Troy,  New  York—  12181 


(lt>  -5b/ y 1  hi 


Sr 


J 


>- 

cu 

o 

C_) 


D  D  C 

nVZ@EIME 


W  FEB  4  1900 


fEtSEUTTE 

B 


.j 


i  L-U 


*This  work  was  supported  by  the  Air  Force  Office 
of  Scientific  Research,  Grant  Number  AFOSR-7 5-2818 . 


DISTRIBUTION  STATEMENT  A 

Approved  for  public  release; 
Distribution  Unlimited 


X  ' 


j 


ABSTRACT 


Collocation  methods  using  both  cubic  polynomials  and  splines 
in  tension  are  developed  for  second  order  linear  singularly- 
perturbed  two-point  boundary  value  problems.  Rules  are  developed 
for  selecting  tension  parameters  and  collocation  points.  The 
methods  converge  outside  of  boundary  layer  regions  without  the 
necessity  of  using  a  fine  discretization.  Numerical  examples 
comparing  the  methods  are  presented. 
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1.  Introduction 


We  consider  the  numerical  solution  by  collocation  methods  of 
the  singularly-perturbed  second-order  linear  boundary  value  problem 
Ly  =  ey''  +  p(x)y'  +  q(x)y  =  f  (x)  ,  a  <_  x  £  b  , 

any(a,e)  +  ct^y'U/e)  *  A,  a21y(b,c)  +  a22y'(b,e)  »  B  , 

where  e  is  a  small  positive  parameter.  The  functions  p,q,  and  f 
are  assumed  to  be  smooth  functions  of  x  on  [a,b]  which,  along  with 
aij'  i'  j  =  1/  2,  and  may  depend  on  e  provided  they  are 
bounded  as  e  -*■  0.  We  further  assume  that  (1.1,2)  has  a  unique 
solution  on  [a,b]  for  all  e  sufficiently  small. 

The  problem  (1.1,2)  has  been  intensively  studied  analytically 
(cf.  Cole  [6],  Eckhaus  [11],  or  O'Malley  [18])  and  it  is  known 
that  its  solution  generally  has  a  multiscale  character,  i.e.,  it 
features  regions  called  "boundary  layers"  where  the  solution 
varies  rapidly.  Away  from  the  boundary  layers,  the  solution  is 
approximately  determined  by  neglecting  the  ey' '  term  in  (1.1) 
and  perhaps  one  or  both  of  the  boundary  conditions  (1.2). 

The  problem  has  also  been  extensively  studied  numerically, 
and  it  is  known  that  most  classical  methods  fail  when  e  is  small 
relative  to  the  mesh  width  h  that  is  used  for  the  discretization 
of  the  operator  L.  There  are,  however,  several  finite  difference 
methods  (cf.  Abrahamson,  Keller,  and  Kreiss  [1]  ,  Berger,  et  al.  [4] , 

II' in  [16],  Kreiss  [17],  and  Pearson  [19,20]),  Galerkin-f inite 
element  methods  (cf.  Hemker  [15],  de  Groen  and  Hemker  [10],  and 
Heinrich,  et  al.  [13,14]),  and  methods  based  on  singular  perturba¬ 
tion  theory  (cf.  Flaherty  and  O'Malley  [12]  and  Steele  [28])  that 
do  not  require  h/e  to  be  small.  A1K  office  of  scientific  research  (afsc) 
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Our  aim  in  this  paper  is  to  show  that  collocation  methods 
with  either  piecewise  polynomials  or  splines  in  tension  can  furnish 
accurate  numerical  approximations  of  (1.1,2)  when  either  h/e  is 
small  or  large.  Splines  in  tension  were  first  used  by  Schweikert  [26] 
as  a  means  of  eliminating  spurious  oscillations  in  curve  fitting 
with  cubic  splines.  They  have  been  subsequently  studied  by  Cline 
[5],  Pruess  [22],  Spath  [27],  and  de  Boor  [9,  Chap.  16].  Between 
knot  points  a  spline  in  tension,  or  a  taut  spline  in  the  language 
of  de  Boor  [9],  is  an  L-  spline  satisfying  the  differential  equation 
(1.3)  (T' 1  -  p2r) ' •  =  0  , 

subject  to  appropriate  continuity  conditions  at  the  knots.  The 
quantity  p  is  called  the  tension  parameter.  When  p  =  0,  t  is  a 
cubic  polynomial  spline;  however,  when  p  >  0,  the  solution  of 
(1.3)  is  a  linear  combination  of  the  four  functions  1,  x,  epx,  e~px. 
The  exponential  functions  should  be  better  suited  than  polynomials 
at  following  the  rapid  variations  that  are  typically  found  in 
singular  perturbation  problems.  Indeed,  exponentials  have  been 
used  in  the  Galerkin  methods  of  Hemker  [15]  and  de  Groen  and  Hemker 
[10],  II' in's  difference  method  [16],  and  the  singular  perturbation 
methods  of  Flaherty  and  O'Malley  [12]  and  Steele  [28]. 

Equation  (1.3)  is  the  same  as  that  for  the  transverse  deflec¬ 
tion  of  a  classical  Euler-Bernoulli  elastic  beam  that  is  subjected 

2 

to  a  tensile  force  proportional  to  p  ;  hence,  the  name  spline 
in  tension. 

In  Section  2  of  this  paper,  we  construct  a  basis  for  the 
tnesion  splines  and  use  it  to  obtain  the  collocation  equations. 

In  Section  3,  we  obtain  asymptotic  approximations  to  the  solution 
of  (1.1,2)  in  the  two  special  cases  when  ]p(x)  |  >_  p  >  0  and  when 
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p(x)  :  0  on  [a,b]  and  use  these  to  select  tension  parameters. 

In  Section  4,  we  discuss  the  selection  of  collocation  points  that 
are  in  some  sense  optimal  and  present  some  formal  error  estimates 
in  regions  not  containing  boundary  layers.  In  Section  5,  we 
apply  our  methods  to  some  examples;  and  in  Section  6,  we  discuss 
the  results. 

Not  surprisingly,  the  results  indicate  that  taut  splines 
provide  better  approximations  within  boundary  layers  and  poly¬ 
nomials  provide  better  approximations  elsewhere.  This  suggests 
the  possibility  of  applying  tension  only  within  boundary  layers. 

This  would  either  require  an  a  priori  knowledge  of  the  location 
of  the  boundary  layers  or  an  automatic  procedure  for  finding  them. 

A  tentative  procedure  for  automatically  locating  boundary  layers 
is  presented  in  Section  5. 

We  anticipate  that  our  methods  would  be  useful  on  other  pro¬ 
blems,  such  as  initial-boundary  value  problems  for  parabolic  partial 
differential  equations  involving  diffusion,  convection  and/or 
reaction. 

2.  Collocation  Equations  and  the  Tension  Spline  Basis 

In  the  usual  method  of  collocation  (cf.  Ascher,  Christiansen, 
and  Russell  [2],  de  Boor  and  Swartz  [8],  Russell  [23],  or  Russell 
and  Shampine  [25]),  one  introduces  a  partition 

1)  As  i  {a  =.  xQ  <  Xl  <  .  .  .  <  xN  -  b) 

of  [a,b]  into  N  subintervals  and  approximates  the  solution  of 
(1.1,2)  by  piecewise  polynomials  y^(x)  £  M(AN,k,m),  where 

2)  M(AN,k,m)  =  {weCm[a,b]  |  wrestr  tQ  T  eP^I^;  i  =  1,2,...,N}  . 

Here 
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and  P^(E)  denotes  the  class  of  polynomials  having  at  most  degree 
k  on  E.  The  dimension  of  M(An,k,m)  is  N(k-m)  +  m  +  1  and  one 
determines  y^(x)  by  collocating  at  N(k-m)  +  m  -  1  points 
z^  e  [a,b],  i.e.,  by  enforcing 

(2.4)  Ly^(z^)  =  f(z^)  ,  i  =  1 , 2 , . . . ,N (k-m)  +  m  -  1  , 

and  by  requiring  y^  to  satisfy  the  boundary  conditions  (1.2). 
Convergence  and  the  order  of  accuracy  of  these  methods  depend  on 
m,k,AN,  and  the  choice  of  and  are  discussed  in,  e.g.,  [8]. 

However,  Ly^  should  exist  and  necessary  continuity  conditions  on 
M(AN,k,m)  are  k  2  and  m  >_  1.  In  addition,  de  Boor  and  Swartz  [8] 
show  that  the  maximal  order  of  convergence  in  the  largest  sub¬ 
interval  length  is  achieved  by  selecting  the  collocation  points 
as  an  appropriate  number  of  Gauss-Legendre  points  on  each  subinterval. 

Unfortunately,  collocation  at  the  Gauss-Legendre  points  with 
piecewise  polynomials  is  knov’n  to  behave  rather  poorly  on  singularly- 
perturbed  problems  for  any  partition  where  e  is  much  smaller  than 
the  minimum  subinterval  length  (cf.  Hemker  [15]  and  our  Example  1 
in  Section  5) .  It  would  be  overly  restrictive  and  in  most  cases 
impractical  to  require  a  partition  with  subinterval  lengths  of 
order  e  and  we  seek  to  avoid  this  situation  by  changing  the  loca¬ 
tions  of  the  collocation  points  and/or  adding  exponential  functions 
to  M(A  ,k,m).  Thus,  we  also  consider  approximations 
vh(x)  e  E(A  , k,m,p),  where 

E(AN,k,m,p)  =  {w  e  cm [a, b] | wrestr _  to  L  e  span  ( ?k ( I ± ) , 
p.x/h.  -o-x/h. 

(2.  5)  e  ,  e  x)  ;  i  =  1,2,  .  .  .  ,N} 


(2.  6) 


h^  =  x^  -  x^_^,  i  -  1,2,. ..,N. 

The  quantities  p^,  i  =  1,2,...,N,  are  called  tension  parameters 
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and  will  subsequently  be  selected  to  approximate  the  rapidly 
varying  part  of  the  exact  solution  of  (1.1,2). 


In  this  paper,  we  have  confined  our  attention  to  collocation 
with  piecewise  cubic  polynomials  belonging  to  M(AN,3,1)  and  splines 
in  tension  or  taut  splines  belonging  to  E(AN,l,l,p),  which  are 
both  spaces  of  dimension  2(N+1).  For  reasons  of  computational 
convenience  and  efficiency,  it  is  usual  to  construct  a  basis  for 
M(A„,3,1)  that  satisfies  the  C^[a,b]  continuity  requirement  and 
where  each  element  has  support  extending  over  only  two  subintervals 
We  proceed  similarly  for  E(AN,l,l,p);  thus,  we  write  y^(x)  in  the 
form 


!2.7) 


yh(x) 


N 

E  (c  .t.  (x,p)  +  d.cr .  (x,p)  ] 
i=0  11  ii~ 


where  the  basis  t^(x,P),  o^(x,P),  i  =  0,1,. . .  ,N,  is  defined  in 
terms  of  the  "canonical  basis  elements"  nQ  and  n 1  as  follows: 


(2.  8 )  t ^  (x,  p)  H 


X  . 

VHq~' '  pi]  '  xetx.^^.l 


x-x . 


V-K^  .  °i+l)  '  X  e  (Xi'Xi+l’  ' 


,  otherwise 


and 


{2.9)ai  (x,  p) 


x .  -x 

“hinl  {“H~  '  pi)  '  x  £  [xi-l'xi] 


x-x 


-  “w1  -  x  £  lxi'xi+i> 


,  otherwise 


The  functions  ru(s,p)  and  n-i(s,p)  are  defined  on  0  <  s  <  1  as 


(2.11) 


n-L  ( s,  p ) 


[2s  -  1 


sinh  p/2  (2s-l) ,  ,  1*  r  -i  cosh 

- sinh  p/2  1  +  2K1  11  “  cosh  p/2  1 


where 

(2.12)  Kq  =  1/pu)  (p/2)  ,  Kx  =  (1/P) coth  P/2,  co(z)  »  coth  z  -  1/z  . 

Observe  that  Hq  and  each  satisfy  the  differential  equation 

(2.13)  (n£'  -  P2nk)"  =  0,  0  <  s  <  1  ,  k  =  0,1, 

subject  to  the  boundary  conditions 

I  t 

(2.14)  nk(o,p)  =  sQk  ,  okd,p)  =  o  ,  nk(0,p)  =  6lk  ,  nkd#p)  =  o  ,  k  =  o, 

I 

where  5  ,  denotes  the  Kronecker  delta  and,  in  this  context,  (  ) 
lk 

denotes  differentiation  with  respect  to  s.  Using  (2.14)  and  (2.7-9), 

we  see  that  c.  =  y.  (x.)  and  d.  =  dy,  (x.)/dx  ,  i  =  0,1,..., N. 
l  n  l  i  n  i 

As  p  tends  to  zero  nQ  and  approach  the  usual  canonical  basis 

elements  for  M(A  . 3,1),  i.e.,  that  of  a  cubic  Hermite  interpolating 
N 

polynomial  on  0  <  s  <  1,  Thus, 

(2.15)  n0(s,p)  =  (1-s)  2  (l+2s)  +  0(p2)  ,  Ti1(s,p)  -  s(l-s)2  +  0(p2)  . 

For  large  values  of  p  rig  and  n-j_  become 

(2.16)  nQ(s,p)  -  1  -  s  -  [  2s-l+e”pS  -e_p  (1's)  ]  +  0(e_p/p) 


Vs' p>  (lW°Sl  - 

Thus,  in  the  interior  of  (0,1)  rig 
by  the  linear  functions 

rig  (s,  p)  -  (1-s)  -  (2s-l)/  (p-2) 


1 

P (P-2) 
and  n 


[s-e“p(1~s) ]  +  0(e-p/p)  . 

^  are  asymptotically  given 

( s ,  p )  ~  ( 1— s )  /  p  -  (2s-l)/o(p-2) 


Both  rig  and  converge  uniformly  on  0  £  s  <_  1  as  p  -*■  »  to  1-s  and  0, 


w 


respectively;  however,  their  derivatives  exhibit  boundary  layer 

behavior  at  s  =  0  and  1.  Since  x^  and  ck  (via.  (2.8,9))  behave 

similarly  at  the  knots  x^^,  x^,  x^+1,  the  numerical  approximation 

yh  may  have  internal  boundary  layer  jumps  of  height  0 ( 1/p )  even 

when  the  exact  solution  is  smooth.  We  demonstrate  this  phenomena 

in  example  1  of  Section  5.  The  functions  Hq  and  are  plotted 

for  a  small  and  a  large  values  of  p  in  Figures  la  and  lb,  respectively. 

A  discrete  system  for  determining  c^,d^,  i  =  0,1,..., N  is 
obtained  for  a  given  set  of  tension  parameters  , i=l, 2, . . . ,N  by 
collocating  at  2N  points  z^, i=l, 2, . . . ,N  on  [a,b]  and  by  satisfying 
the  boundary  conditions  (1.2).  For  simplicity  we  place  two  collo¬ 
cation  points  symmetrically  disposed  on  each  subinterval,  i.e., 

(2.17)  Z-.  ,  =  x.  .  +  t.h.  ,  z_.  =  x.  .  +  (l-t.)h.  ,  i  =  1,2,. ..,N, 

2i-l  i-l  li  2i  i-l  ii 

for  an  appropriate  choice  of  ^  e  [0,1/2).  Then,  using  (2.17), 

(2.7-9),  and  (1.1)  in  (2.4)  we  find  the  discrete  system  on  1^  to 


be 


(2. 18) 


'il 


(ti> 

Ai2(ti> 

*i3(ti> 

£i4 (ti} 

(1-t.) 

ii2(1“ti) 

*13 {1"ti) 

ti4(1-ti) 

where 


c .  . 
i-l 


li-l 


c  . 

i 


d. 

i 


i  =  1 , 2 , . . . ,N 


f . (t.) 

l  l 


f, (1-t.) 


L  1 


1  J 


(2.19) 


(t)  =  (c/hi  )nQ  (t,pi)  +  (pi  (t)/hi)  n0  (t,  pi)  +  qi  ( t)  n0  (t,  pi)  , 

^i2  (t )  =  (e/hi)ri1  (t,pi)  +  p.^  { t)  n1  ( t,  pi>  +  qi  (t)  hin1  ( t,  pi)  , 
ii3(t)  =  (e/hi2)nQ  (l-t,pi)  -  (pi  (t) /hi)  nQ  (1-t,  pi)  +  ( t)  nQ  (1-t,  p±) 


*i4(t)  -  -[(e/hi)n1  (l-t,pi)  -  p^.  (t)  (1-t,  pi)  +  qi  ( t)  r)1  ( 1-t,  p^  ] 
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(2.20) 


(2.21) 


and  f^(t)  =  f(x^_^+th),  etc.  Substituting  (2.7-9)  into  (1.2) 
gives  the  boundary  conditions 

allc0  +  a12d0  =  A  '  a21cN  +  **22^  =  B  * 

Thus,  the  2(N+1)  dimensional  discrete  system  (2.18-20)  has  the 
following  structural  form 


X 

X 

V 

r  A 

X 

X 

X 

X 

d 

o 

X 

X 

X 

X 

C1 

X 

X  X 

X 

dl 

• 

X 

X  X 

X 

• 

2S 

• 

•  • 

X 

X 

X 

X 

• 

• 

yv 

X 

X 

X 

X 

CN 

fN(1~tN) 

X 

X 

dM 

B 

_ 

_ 

L  J 

/ 


where  each  x  denotes  a  nonzero  entry.  We  solve  (2.18-20)  for 
prescribed  values  of  p.^  and  t^,  i  =  1,2,...,N,  by  an  alternating 
row  and  column  pivoting  algorithm  due  to  Varah  [29].  This  procedure 
is  numerically  stable  and  requires  no  storage  additional  to  that 
needed  for  the  nonzero  entries  in  (2.21). 

3.  Asymptotic  Solutions,  Green's  Functions,  and  the  Selection  of 
Tension  Parameters 

In  this  section,  we  present  asymptotic  approximations  of  the 
solution  of  (1.1,2)  and  of  its  Green's  function.  They  will  be 
used  to  select  the  tension  parameters  and  in  Section  4  to  select 
collocation  points.  We  shall  not  attempt  to  do  this  in  all  gener¬ 
ality,  but  rather  by  considering  two  special  cases  of  the  problem 


8 


i  i  i 

(3.1)  Ly  =  ey  +  p(x)y  +  q(x)y  =  f  (x)  ,  a  ^  x  <_  b  , 

(3.2)  y(a,e)  =  A  ,  y(b,e)  =  B  , 

when  (Problem  1)  j.p(x)  |  >_  p  >  0  and  (Problem  2)  when  p(x)  =0, 
q(x)  <  q  <  0  for  x£[a,b].  In  either  case,  any  boundary  layers 
are  at  the  ends  of  [a,b];  thus,  there  are  no  turning  points  and 
no  interior  nonuniformities.  We  consider  such  problems  among  the 
examples  of  Section  5. 

The  Green's  function  G(x,£)  associated  with  the  operator  L 

i 

and  homogeneous  boundary  conditions  (3.2)  on  the  interval  [a,b] 
satisfies 


(3.3) 

L*G  (x, £ )  =  gG^^  -  (p(5)G)? 

+ 

q(5)G  =  0, (a,x)  u  (x,b) , 

(3.4a) 

G  (x,  a)  =  G  (x,  b)  =  0 

/ 

(3.4b) 

G(x,x+)  -  G(x,x  )  = 

0 

t 

(3.4c) 

G^  (x,x+)  -  G^  (x,x“) 

= 

1/e  , 

where  the  subscript  5  denotes  partial  differentiation. 

We  use  the  WKB  method  to  construct  our  asymptotic  approximations 
of  the  solutions  of  (3.1,2)  and  (3.3,4).  Since  the  details  of 
this  method  are  well  known  (cf.  Wassow  [30]),  we  only  present  the 
results  and  omit  their  development. 

3.1.  Problem  1:  lp(x) 1  >  p  >  0  for  xe  [a,b] 

We  consider  the  case  when  p(x)  >  0  on  [a,b] .  The  case  when 
p(x)  <  0  is  handled  in  an  analogous  manner.  Using  the  WKB  method, 
we  find  the  following  0(e)  approximations  to  the  two  fundamental 
solutions  of  (3.1)  (cf.  Hemker  [15]): 

3.  5a,  b)  Y  (x,  \  )  =  exp^/^q  ( z)  /p  (z)  dz }  ,  11(5, x)  =  exp{-/x  (p  ( z) /e-q  ( z) /p  ( z)  )  dz  . 

The  solution  of  (3.1)  satisfying  the  boundary  conditions  (3.2)  is 
given  by 
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(3.6a)  y(x,e)  =  YR (x )  +  [A  -  YR(a)]II(a,x)  +  0(e) 

where 

(3.6b)  Y_  (x )  =  BY  (x,  b)  -  /b(f  (z)/p(z)  )Y(x,z)dz  . 

R  x 

The  term  [A  -  Y_  (a)  ]  II  (a,x)  is  exponentially  small  outside  of  a 
boundary  layer  of  width  0(e)  near  x  =  a.  For  a  <  x  <_  b  y(x,e)  ~  YR(x) 
where  Y_(x)  is  the  solution  of  the  reduced  problem 

K 

(3.7a)  p(x)YR(x)  +  q(x)YR(x)  =  f  (x)  ,  a  <_  x  <_  b  , 

(3.7b)  YR(b)  =  B  , 


(3.8a) 


I  t 

obtained  by  neglecting  the  ey  term  in  (3.1)  and  the  boundary  con¬ 
dition  (3.2)  at  x  =  a.  The  problem  with  p(x)  <  0  has  a  solution 
with  a  boundary  layer  near  x  =  b  and  a  reduced  solution  satisfying 
(3.7a)  subject  to  the  initial  condition  YR(a)  =  A. 

In  a  similar  manner,  an  0(e)  approximation  to  G(x,£)  satisfying 
(3.3,4)  is  found  as 

G(x,£)  =  a(x)  (II(£,b)Y*(a,b)  (II(a,x)  -  Y*(x,a)] 

f  II(£,x)  ,  a  <  5  <  x 


-H  (a,x)  Y*  (a,  £)  +  1 


}  +  0(e) 


Y*(X,S)  ,  x  <  £  <  b 


where 

(3.8b)  a(x)  =  -  p(x)/[p2(x)  -  2eq(x)  +  sp  (x)  ]  , 

(3.8c)  Y*  (x,  £ )  =  exp  {  /^  ( q(z)  -  p  (z))/p(z)dz  , 

x 


and  n  (£,x)  is  as  in  (3.5b).  As  a  function  of  £,  G(x,£)  has  boundary 
layers  at  £  =  b  and  £  =  x  . 

3.2.  Problem  2:  p(x)  =  0  ,  g(x)  <_  q  <  0  for  xe[a,b]  . 

In  this  case,  the  WKB  method  gives  the  following  0(vT)  approxi¬ 
mations  to  the  two  fundamental  solutions  of  (3.1)  (cf.  Hemker  [15]). 

( 3  .  9a ,  b)  H1(x,£)  =  q  (x)  (x,  £)  ,  H2(x,£)  =  q  (x)  "1/4H  ( £  ,  x) 
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where  for  this  problem 


(3.9c)  n  (£,x)  =  exp  {-/X  /-q(z)/e  dz}  . 

£ 

The  character  of  the  solution  depends  critically  on  the  sign  of 
q(x).  When  q(x)  <  0,  the  solution  is  exponential;  and  when  q(x)  >  0, 
the  solution  oscillates  rapidly  with  period  2ir/e/q (x)  .  We  would 
not  expect  taut  spline  approximations  to  be  useful  in  the  oscillatory 
case  and,  thus,  we  confine  our  attention  to  problems  with  q(x)  <  0 
on  [a,b].  The  use  of  "splines  under  compression”  for  oscillatory 
problems  is  currently  under  investigation  by  Coyle  and  Flaherty  [7] . 

Using  (3.9),  the  solution  of  (3.1,2)  is  given  as 
(3.10)  y(x,e)  =  [A  -  f  (a)/q(a)  3  (q(a)/q(x)  ]  1/,4II  (a,x) 

+  [B  -  f  (b)/q(b)]  (q(b)/q(x)]1/4n(x,b)  +  f(x)/q(x)  +  0(/e)  • 
Outside  of  the  boundary  layers,  which  extend  over  0(/E)  neighbor¬ 
hoods  of  x  =  a  and  b,  the  solution  y(x,e)  -  YR(x),  where  YR(x)  is 
the  solution  of  the  reduced  problem 

q (x) Yr (x)  =  f  (x)  , 

1  I 

obtained  by  neglecting  the  ey  term  in  (3.1)  and  both  boundary 
conditions  (3.2). 

Problem  2  is  self-adjoint,  so  the  WKB  approximations  (3.9) 
can  also  be  used  to  construct  the  following  0 (1)  approximation  to 
G(x,£)  satisfying  (3.3,4): 


(3.11) 

g(x,o 

=  |[e2q(x)q(0  ]"1/4(n( 

a ,  x)  II  (a ,  £ )  +  n(x,b)H(£,b) 

j"  n U,x)  ,  a  £  C  1  x 

•—  « 

1 

+  0(/e)  }  . 

L  n  (x,  £)  ,  x£5<_b 

As  a  function  of  5,  G(x,£)  has  boundary  layers  on  both  sides  of 
£  =  x,  and  is  unbounded  as  0(l//e)  as  e  -*■  0. 


3.3.  Selection  of  Tension  Parameters 


We  want  the  tension  parameters  to  approximate  the  rapidly 
decaying  solutions  (3.5b)  or  (3.9)  of  Problems  1  or  2,  respectively/ 
and  so  we  choose  on  the  subinterval  1^  as 


(3.12a)  p 


|p(xk)/e  -  q(xk)/p(xk)  | 


if  jp(xi_1)+p(xi) |/e  >  |q(xi_1)+q(xi)  | 

# 


/-  [q(xi_1)+q(xi)  ]/2e 


if  |p(xi_1)+p(xi) |/e  <  |q(xi_1)+q(xi) 


i  -  1/ 2, . . . ,N, 


where 


(3.12b)  k 


i  -  1  ,  if  Cp(xi_1)  +  p(xi)]/s  >  0 


[  i  /if  fP(xi_1)  +  p(xi)]/e  <  0 

Thus,  for  Problem  1  with  p(x)  >  0  on  [a,b] 

(3.13a)  Pi  =  hi!p(xi_1)/£  -  q(xi_1)/p(xi_1)  |  , 

and  for  Problem  2 

(3.13b)  p±  =  h±  /-[q(xi_1)  +  a\xi)]/2e  . 

However,  we  use  (3.12)  computationally  even  when  the  conditions 
of  Problem  1  or  2  are  not  satisfied,  e.g.,  when  there  are  turning 
points. 

The  solution  of  the  collocation  equations  (2.18-20)  with  the 
tension  parameters  specified  by  (3.12)  will  give  the  exact  solution 
of  (1.1,2),  for  any  choice  of  t^  e [0,1/2),  whenever  f (x)  is  a  linear 
polynomial  and  either  1)  p(x)  =  p  and  q(x)  =  0  or  2)  p(x)  =  0  and 
q (x )  =  q  <  0.  This  is  because  the  solutions  of  these  problems 
are  elements  of  the  approximating  space  E(AN,l,l,p). 

We  close  this  section  by  applying  the  method  of  collocation 
with  splines  under  tension  to  the  example 
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(3.14)  Ly  =  ey  +  p(x)y  =  f  (x)  ,  a  <_  x  <_  b,  y(a)  =  A,  y(b)  =  B, 

with  p(x)  >0  on  [a,b].  Using  (3.5,6)  the  solution  of  this  problem 


(3.15a) 


y(x,e)  =*  YR(x)  +  (a  -  YR (a)  ]  exp  {-/( p(z)/e)dz} 


where 


(3.15b)  Y  (x)  =  B  -  /“(f  (z)  /p(z)  )dz  . 

R  x 

For  the  present,  we  choose  t^  =  0  and  then  use  (2.10-12)  and  (2.19) 
in  (2.18)  to  obtain  the  discrete  system 

(epi/hi)  (-^i  -  |ydj  u)"1(pi/2)  +  |vd±  coth  Pi/2  +  p(xi_1)di_1  =  f(xi_1) 
1.16)  (epi/hi)  (^i  -  |udi)  ai-1(pi/2)  -  |vdi  coth  p±/2  +  p(x.)d.  =  Ux±)  , 


i  —  1,2,. ..,N 


°0  *  &  '  CH  =  B 


where 


(3.17)  V(  ).  =  (  )•  -  (  )w  ,  U(  )•  H  (  ).  +  (  )._;L  , 

and  ui(z)  is  defined  by  (2.12).  Using  (3.13a)  we  select  p^  =  h^p(x^_1)/e 
and  assume  that  the  partition  has  been  chosen  so  that  pi  >>  1, 
i  =  1,2,...,N.  In  fact,  suppose  that  p^  is  large  enough  tc 
approximate  oj(p^/2)  and  coth  p^/2  by  (1-2/p^)  and  1,  respectively. 

Then  (3.16)  become 

(e/hi)  (2Vci/hi  -  ud^  +  P  (xi-1)  Vci/hi  =  f  (xi_1)  , 

(3.18) 

(pp(xi))di  =  u  f (x i )  ,  i  =  1,2,...,N,  cQ  =  A,  cN  =  B 

Thus,  the  solution  is  approxiamtely  determined  as  the  solution  of 
3.19a, b)  cN  =  B  ,  p(xi_1)  7ci/hi  =  f  (x^)  +  OU/hu)  ,  i  =  N,  N-l,...,2, 

(3.19c)  (yp(xi))di  =  uf (xi)  ,  i  =  N,N-1, _ ,1, 


13 


19d,e)  cQ  =  A  ,  dg  =  -(p(a)/e)  [cQ  -  (a) /p  (a)  +  0(e/h1)]. 

Equations  (3.19a,b,c)  can  be  recognized  as  0(h)  (where  h  is  the 
maximum  subinterval  length)  "upwind"  difference  approximations  to 
the  reduced  problem,  while  (3.19e)  gives  the  initial  slope  of  the 
solution  in  the  boundary  layer  correct  to  0(h/e). 

4.  Selection  of  Collocation  Points 

Our  aim  in  this  section  is  to  suggest  some  special  choices 
of  collocation  points  that  may  be  used  to  reduce  the  errors  in  methods 
for  Problems  1  and  2.  We  confine  our  attention  to  the  two  limiting 
cases  of  zero  tension  (p^  =  0)  and  large  tension  (p^  >>  1)  for 
i  =  1,2,...,N.  For  simplicity,  we  consider  uniform  partitions 
with  h  =  hu,  i  =  1,2,...,N,  and  assume  that  hp(x)/e  >>  1  for 
Problem  1  and  h/-q (x) /e  >>  1  for  Problem  2.  No  detailed  rigorous 
error  analysis  will  be  given;  however,  some  formal  error  estimates 
are  obtained  on  subintervals  not  containing  boundary  layers. 

4.1.  Error  Formulas 

It  is  well  known  (cf.  [8],  [23],  or  [24])  that  the  pointwise 
error  in  collocation  methods  for  (1.1,2)  satisfies  an  equation  of 
the  form 

(4.1)  e(k)  (x)  =  y(k)  (x)  -  y^k)  (x)  =  fb  G(x,£)r(£)d5  ,  k  =  0,1 

n  a  3xk 

where  the  residual 

(4.2)  rU)  =  Ly(£)  -  Lyh(S)  =  f(C)  -  Lyh(C)  • 

It  is  convenient  to  introduce  the  local  transformation 

(4.3)  £  =  +  hs  ,  0  <  s  <  1 

/v 

on  the  subinterval  1^  and,  as  in  Section  2,  let  f^ts)  =  f (xi_1+hs) , 
etc.  Then  (4.1)  becomes 
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(4’.  4a) 


(4.4b) 

(4.5a) 

(4.5b) 

(4.6) 


(4.7) 

(4.8) 


(4.  9) 

(4.10) 


00 


N 


(x)  *  h£  Z1 


i=l  0  3x 


— £  G (x,xi_1+hs)  ri(s)ds,  k  =  0,1, 


where 


ri<s)  -  f±(s)  -  Liyh_(s)  , 


A  A 


Liyi(s)  =  (£/h2^i  {s)  +  LR.yi(s)  ' 


and  Lr  is  the  reduced  operator 

i 


A  A 


A  I 


LR  yi(s)  "  (Pi (s)/h) y^_ (s)  +  q^sjy^s) 


Let  Pf^(s)  be  the  linear  interpolant 

Pfi(s)  =  [(l-ti-s)fi(ti)  +  (s-ti)fi(l-ti)]/(l-2ti) 

A 

to  f.(s)  at  the  two  collocation  points  s  =  t.  and  1  -  t.  on  I . . 

1  c  1X1 

A  A  A 

Since  the  collocation  equations  (2.4,17)  imply  L^y^  *  f .  at  s  8  t., 

A  A  A  1 

1  -  t.,  we  have  PL.y.  =  Pf .  and  (4.4a)  may  be  written  as 
X  1  1 

(\C\  ^  1  aK  A 

e'  '  (x)  *  h£  Z  — 5-  G  (x,x .  ,+hs)  (l-P)r.  (s)ds  ,  k  =  0,1  . 
i=l  0  3xK  1“1'  1 

The  interpolation  error 

(1-P)ri(s)  *  (s-t^  (s-l-ti)ri(ti,  1-t^  s]  , 


where  r^ [s^, s^, . . . , s^l  denotes  the  k  th  divided  difference  of  r^ 

A 

at  the  points  sQ, s^, . . . , s^.  This  form  of  (l-P)r^(s)  suffices  when 
Pi  =  0;  however,  when  >>  1  a  more  detailed  form  is  needed.  In 
this  case,  we  assume  that  p .  is  large  enough  to  neglect  terms  of 

-Pf/2  1 

0(e  )  relative  to  unity  and  use  the  large  tension  approximations 


(2.16)  and  (2.7-9)  in  (4.5)  to  get 

A.A.  P.  “P'S  A  —  P  .(l-S) 

Liyh ,  *  —  te  »i  “i(s>  '  e 


YiV.tsn^.Yh.ls)  , 


where 


8i  =  Vc./h  -  di_1  -  (Vdi)/pi  ,  Yi  =  Vc./h  -  &L  +  (Vdi)/pi 
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(4.11) 


(4.12) 


(4.13) 


(4.14) 


ui(s)  =  ep.j/h  -  Pi(s)  +  qi(s)h/pi,  vi(s)  »  ep^h  +  pi(s)  +  qi(s)h/pi  , 

and  Y.  (s)  is  the  linear  polynomial  part  of  y.  (s) ;  thus, 
fti  ni 

Vhi(s)  *  p^5  (s)  (Vc^/h  -  udi/pi) 

+  +  di_1h/Pi)  (1-s-l/Pj^)  +  (ct  -  d1h/pi)  (s-l/Pj)  ]  . 

A  A 

The  choice  of  p^  given  by  (3.12)  makes  u^(s)  =  0(h)  when  p^(s)  >  0 

A  /S  ^  /\  "  2 

on  1^,  v^(s)  =  0(h)  when  p^(s)  <  0  on  1^,  and  u^s)  =  v^(s)  *  0(h  /p^) 

A 

when  p^(s)  =  0  on  1^  . 

Using(4.9)  and  (4.4b),  we  have 


A  A 


(l-P)r.(s)  =  (s-t.)  (s-l+t.)  {f  .  [t.,l-t.,sj  -  L_  Y.  [  t .  ,  1— t . ,  s  ]  } 
l  i  1111  h^  l  i 


--j  (6iui(s)g(s,ti)  -  Yivi(s)g(l-s,ti) 


+  ( s— t  ^ ) (s-l+ti) 


-p . t .  -p  . (1-t . ) 
ii  l  i 

e  -e 

l-2t . 
l 


((.■.ll-tj.sItVill-ti.sl) 


I  - 

-(s-ti) ( s— 1+t^) y e  Siui[ti,l-ti,s] 


-Pi (1-t.)  ~ 

+e  Yivi(ti,l-ti,s] 


where 

-PdS  -P4td  -p. (1-t.) 

g(s,ti)  =  [(l-2ti)e  1  +  (s-l+ti)e  1  1  -  (s-t^e  1  1  ]/(l-2t.) 

-0./2 

The  assumption  that  terms  of  0(e  )  are  negligible  will  specifically 

allow  us  to  drop  all  terms  in  (4.13,14)  involving  the  factor 
-Pi  d-t L) 

e  since  t^  <  1/2.  This  will  be  done  in  all  further  uses 

of  (4.13,14)  except  where  noted. 

A 

It  remains  to  use  the  formulas  (4.8)  or  (4.13)  for  (l-P)r^(s) 
together  with  the  approximations  (3.8)  or  (3.11)  for  the  Green's 
functions  in  (4.7)  and  find  appropriate  choices  for  collocation 
points.  One  problem  is  that  the  errors  given  by  (4.7)  depend  on 
the  unknown  numerical  solution  y^.  This  would  not  be  a  serious 
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difficulty  if  yh  and  y^  were  bounded  as  e  +  0  for  fixed  h.  We  have 

I 

shown  by  example  in  Section  3  that  c^  and  d^  (hence,  y^  and  y^) 
are  bounded  away  from  the  boundary  layer  region  for  e/h  <<  1  when 
t^  »  0  and  is  selected  according  to  (3.12).  It  is  reasonable 
to  assume  that  this  remains  so  when  t^  is  sufficiently  small;  however, 
it  is  also  relatively  easy  to  show  that  d^  can  be  unbounded  at 
every  knot  point  as  e/h  -*■  0  when  p^  =  0  and  t^  are  the  Gauss- 
Legendre  points.  Little  is  known  about  the  behavior  of  the  numerical 
solution  for  other  choices  of  t^  and  p^.  In  this  paper,  we  shall 
not  attempt  to  find  conditions  for  y^  to  be  bounded  as  e  +  0,  but 
rather  we  shall  make  some  suggestions  for  collocation  points  that 
should  generally  reduce  the  error  in  methods  for  Problems  1  and  2. 

f 

We  note  in  passing  that  if  yh  and  y^  were  bounded,  arguments  similar 
to  those  used  by  Pruess  [21]  or  Russell  and  Christiansen  T24]  on 
related  non  singularly-perturbed  problems  could  be  used  to  remove 
the  dependence  on  y^  from  the  leading  order  terms  in  e(x) . 

4.2.  Collocation  Points  for  Problem  1 

We  again  consider  the  case  when  p(x)  ^  p  >  0  for  xe  [a,b]. 

Let  x  be  a  knot  point,  say  x^,  so  that  there  are  no  discontinuities 
in  derivatives  of  the  Green's  function  on  any  subinterval  and  apply 
the  transformation  (4.3)  to  (3.5b)  and  (3.8c)  to  get 

A  A 

(4.15a)  n(xi-l  +  hs'xj)  “  H  (xi,Xj)TTi(s)  =  6ij1ri(s)  ,  i  <  j 

(4.15b)  Y*(Xj,xi_1  +  hs)  =  Y* (Xj,x^_^) £^(s)  ,  i  >  j 

where 

(4.15c)  TT.(s)  =  exp  {-h/^  [p .  (z) /e  -  q  .  (z) /p  .  (z)  ]  dz  , 

s  11 

(4. 15d)  (s)  =  exp  { h/ s [ ( q . ( z )  -  p . (z) ) /p . ( z) ] dz  . 

0  1  11 
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The  last  form  of  (4.15a)  follows  from  our  assumption  that  uerms  of 

-Pi/2 

0(e  )  are  negligible,  hence,  the  boundary  layer  in  II(x.,x.) 


i  3 


at  Xj  is  well  within  subinterval  I y 

Using  (3.8a)  and  (4.15)  in  (4.7)  we  have 


N  ~ 


e(x..)  =  h[Y*(a,b)6j0  -  Y*(xj,b)]SN  -  h6j0  Y*(a,xi_1)Ri 

N  ~ 

(4.16a)  +  ha  (x  . )  S  .  +  h  £  Y*(x.,x.  . )  R .  , 

3  3  i=i+l  3  1-1  1 


(4.  16b) 


e  (Xj)  =  -  h{Y*  (a,b)  [p*/h  -  a  '  (x.. )  /a  (x  ^  )  ]  6  j  Q  +  Y*(x^,b)  SN 

*  ,  N 

+  h  [p./h  -  a  (x  .)/a(x  .)]<S  .n  £  Y*(a,x.  .)R. 

J  J  J  Ju  j_=]_  1 

*  f  N 

i  -  ha  (x  . )  [p./h  -  a  (x . )  /a  (x  . )  ]S  .  +  h  £  Y  (x.,x.  ,)R. 

3  3  3  3  3  i=j+l  x  ]  i  x  i 


where 

t.17,18)  p*  =  h  [p  (Xj  )  /  e  -  q(Xj)/p(Xj)  ]  ,  Y*(x^,xi)  =  a  (x^ )  Y*  (x  ^  ,  x±)  , 


1. 19a,  b)  S.  =  /1tt  (s)  (l-P)r.  (s)ds  ,  R.  =  ,  (s)  (1-P)  r .  (s)  ds  . 

1  0  1  1  1  0  1 

*  * 

We  call  Pj  the  adjoint  tension  parameter  and  note  that  p^  =  pj+l 

when  the  tension  parameters  are  selected  according  to  (3.13a). 

★ 

Observe  that  P^  is  large  whenever  hp(x)/e  is  large  and  this  can 

I 

induce  large  errors  in  e  (x^)  (cf.  (4.16b)).  These  large  errors 
can  be  confined  to  the  boundary  layer  near  x  =  a  if  is  sufficiently 

f 

small.  However,  in  order  for  e  (x^)  to  be  small  within  the  boundary 
layer  R^,  i  =  1,2,...,N,  must  also  be  sufficiently  small.  In  this 
paper,  we  have  concentrated  on  producing  good  approximations  outside 
of  boundary  layer  subintervals. 

a 

We  first  consider  the  taut  spline  approximation  where  (l-P)r^(s) 
is  given  by  (4.13,14).  Expanding  f (x) ,  p(x) ,  and  q(x)  in  Taylor's 


series  about  a  suitable  point  on  1^  would  reveal  that  f ^ [t^,l-t^, s] 
and  ui[ti,  1-t^  s]  are  0(hz )  while  u^l-t^s]  and  v^l-t^s]  are 
0(h).  As  previously  noted,  the  choice  of  p.  given  by  (3.13a)  makes 

/v 

ui(s)  of  0(h)  and 


(4. 20) 


v^  ( s)  =  pi(0)  +  pi(s)  +  0(h/pi)  . 


If  we  further  assume  that  c.,  d./p.,  and  Vc./h  are  bounded  then 

i  i  i  i 

A  A  2 

L_  Y.  [t.,l-t.,s]  is  0(h  ),  and  to  leading  order  (4.13,14)  become 

(4.21)  (l-PJr^s)  ~  Yivi(s)g(l-s,ti) 

This  is  not  surprising  since  this  term  is  due  to  the  presence  of 
-Pi  (1-s) 

the  e  functions  in  the  taut  spline  approximation  and  these 

functions  are  not  present  in  the  exact  solution  of  Problem  1. 


. 22a, b) 


Substituting  (4.21)  into  (4.19)  gives 

^  A  A 

S.  1  Y.  /  IT.  (s)g(l-s,t.)v.  (s)ds  , 

±  0 


Ri  :  Yi  /  Ci(s)9(l-s,ti)v1 


S±  may  be  further  approximated  by  using  Laplace's  method  (cf.  Bender 

/V 

and  Orszag  [3,  Chap.  6]).  The  essential  idea  is  that  tt^s)  is 
exponentially  small  outside  of  a  small  neighborhood  of  s  -  1  when 
hp^(s)/e  is  large;  hence,  the  integrand  in  (4.15c)  may  be  replaced 
by  its  value  at  s  =  1.  *Using  (4.17)  this  gives 

t”Pj(1—s)  ~ 

[4.23)  S.  Z  y  .  /  e  g  (1-s,  t.  )v.  (s)ds. 


Now,  we  see  that  S.  may  be  reduced  in  magnitude  by  selecting  t. 

-p  (1-s}  1 


such  that  e 

we  require* 

1  "Pi  (1-s) 
(4.24)  /  e  1 

0 


is  orthogonal  to  g(l-s,t^),  i.e.,  using  (4.14) 


1  1  l-t.-l/p.  -p.t. 

gU-s^ds  “  Jf  [I7gT7p|  -  1-2^  e  1=1 


If  terms  of  0(l/p^)  are  neglected,  this  implies 
(4.25)  ti  =  (l/pi)  In  (l+p./p*)  ,  i  =  1,2,...,N 


(s)  ds 
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(4.26) 


(4.27) 


(4.28) 


We  refer  to  this  choice  of  t^  as  Method  1.  It  can  only  be  used 

*  * 
when  p^  and  are  such  that  t^  <  1/2.  When  p^  and  p^  are  large 

t ^  =  0(l/p^)  and  collocation  is  performed  near  the  ends  of  each 

subinterval.  Using  (4.14),  (4.20),  and  (4.25)  in  (4.22b)  implies 

Ri  :  Yi?i(l)v.(l)e  Pl  1  I  YiPi(l)?i(l)  . 

If  both  c^  and  d^  were  bounded  outside  of  the  boundary  layer, 
then  Y-  would  be  0(h)  (cf.  (4.10)).  Thus,  from  (4.16)  the  best 

f 

that  we  could  expect  from  Method  1  is  for  e(x^)  and  e  (x^)  to  be 
0 (h) .  The  computational  evidence  in  Section  5  indicates  that 
this  is  the  case. 

A  second  possibility  is  to  select  t^  so  that  given  by 

(4.22b)  is  reduced  in  magnitude.  This  can  be  done  by  requiring 

/1g(l-s,t. )ds  =  0. 

0  1 

Using  (4.14)  this  leads  to 

fci  =  I  [1  '  p~  cosh"’1  sinh  » 

and  we  refer  to  this  choice  of  t.  as  Method  2.  We  retained  the 
-p.(l-t.)  1 

0(e  1  1  )  term  in  (4.14)  when  obtaining  (4.27)  in  order  that 

ti  approach  the  Gauss-Legendre  point  (cf.  (4.33))  as  p^  ■*  0.  When 

Pi  is  large  t^  ~  (1/p^)  In  (p^/2)  and  in  this  case  (3.12),  (4.18), 

and  (4.23)  imply  that 

S .  2  y • P- (1)/P • 
l  '  1*1  '  i 

The  computational  evidence  in  Section  5  indicates  that  S.^  is  not 

small  enough  to  insure  an  accurate  approximation  of  e  (x.)  at  any 

2 

knot  point.  In  fact,  the  indications  are  that  e(x.. )  -  0  (h  )  while 
'  ~  2 

e  (Xj)  -  0 (h  /t)  when  p^  is  large.  Method  2  may  still  be  used  in 
this  case  if  one  is  not  interested  in  predicting  the  slope  of  the 
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W 


solution  as  long  as  is  not  so  large  that  it  causes  the  discrete 
system  (2.18-20)  to  be  ill-conditioned.  When  is  small  t^ 

I 

approaches  the  Gauss- Legendre  point  and  e(Xj)  and  e  (x^)  will 
4 

approach  0(h  )  (cf.  deBoor  and  Swartz  [8]). 

For  polynomial  approximations,  we  use  (4.8)  in  (4.19)  to 

get 


(4.29a) 


S.  1  f  e 
1  0 


-P,  (1-s) 


(s-t^  (s-l+ti)  ri  [ti#  l-tif  s]ds 


(4.29b) 


^  «  /1;i(s)  (s-t.^  (s-l+ti)ri[ti,l-ti,s)ds  , 


where,  once  again,  Laplace's  method  was  used  to  approximate  the 
singular  integral  S^. 


Either  S.  or  R. 


. .  __  _t.  may  be  reduced  in  magnitude  by  selecting  t. 

1  1  -p* (1-s)  1 

such  that  either  e  1  or  1  is  orthogonal  to  (s-t^) (s-l+t^) 


respectively,  i.e.,  by  requiring  either 

■.  —  p  •  ( l~s)  . 

1. 30a, b)  f  e  (s-t . )  (s-l+t.  )ds  =  0  or  /  ( s-t  )( s-l+t  ) ds  =  0 

0  ii  0  1 

The  option  (4.30a)  gives 


(4.31) 


-  u>(p./2) 

L  pi  l+/l-4w(p TTITTp? 


1  "  ■  1 

where  oj(z)  was  defined  in  (2.12).  This  method  is  referred  to  as 
Method  3.  Once  again,  t^  approaches  the  Gauss-Legendre  pcint  as 


0  and  t^  2  1/p^  when  p^  is  large.  Assuming  that  y^  and 


y,  are  bounded  outside  of  the  boundary  layer  region  and  using 
n 

(4.4b)  and  Taylor's  series  expansions  of  f (x)  and  Ly^(x)  about 

x.  in  (4.29b)  leads  to 
i  2 

(4.32)  Ri  :  "  12  ~  Lyh(xi)]"  . 
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(4.33) 


(4.  34a) 

(4.34b) 

4. 34c, d) 


This  in  turn  via  (4.16)  would  imply  that  e(Xj)  and  e  (Xj)  are 
2 

0(h)  at  knot  points  outside  of  the  boundary  layer.  However,  it 

II  II 

is  typically  possible  to  replace  (Lyh(Xj))  by  (Ly(x^))  +  0(h) 

(cf.  Pruess  [21]  or  Russell  and  Christiansen  [24]).  If  this  were 

1  3 

so,  this  R^,  e(Xj),  and  e  (x^)  would  all  be  0(h  )  at  knots  away 
from  the  boundary  layer.  This  is  in  accord  with  the  computational 
results  of  Section  5. 

The  final  possibility  is  to  select  t^  so  that  (4.30b)  is  zero 
and  this  gives  t^  as  the  Gauss-Legendre  point  on  the  interval 
[0,1],  l. e. , 

t±  =  (1  -  l//3)/2  . 

This  choice  of  t^  is  known  to  give  poor  results  when  hp(x)/e  >>  1; 
however,  in  Section  5,  we  show  that  it  may  be  used  outside  of  sub- 

I 

intervals  containing  boundary  layers  provided  that  yh  and  yh  are 
computed  accurately  enough  at  the  ends  of  subintervals  containing 
boundary  layers. 

4.3.  Collocation  Points  for  Problem  2 

Again,  let  x  be  a  knot  point,  say  x^,  and  use  (3.9)  and  (4.3) 
to  write 

A  A 

n(Xi_i+hs,Xj)  =  n(xi,Xj)iri(s)  =  5ij7t(s)  ,  i  <_  j 

A  A 

IHx^x^hs)  =  H  (Xj  ,xi)  Xi  (s)  =  6ijXi(s)  ,  i  >  j  , 
where  now 

^,(s)  =  exp  {-h/1/-^, (s) /eds) ,  X (s)  =  exp  {-h/s/-§ .  . (s) /e  ds  . 
s  0  1 

The  last  forms  of  (4.34)  again  follow  from  our  assumption  that 
-Pi/2 

terms  of  0(e  )  are  negligible.  Using  (3.11)  and  (4.34)  in 

(4.7)  we  have 
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:4.35a) 


(4. 35b) 

. 36a,b) 

(4.37a) 

(4.37b) 


(4.  38) 


e(x.)  =  j  («oj  S,  +  5Nj  SN  -  S.  -  Sj+1)  , 

e  =  2  +  W  50j  S1  +  (Pt/h  "  qN//qN)  5Nj  SN 

+  ( P*/h  +  q j/ q j )  sj  "  (Pj/h  -  qj/qj  >  ;)+l>  ' 


where 


qj  =  q(Xj) 


Pj  =  h/-q  (x  j  )  / e 


and 


S  .  =  /  it  .  (s)  (e  q  (x . )  q .  (s)  ] 
J  0  J  J  J 


-1/4 


(1-P) r j ( s) ds  , 


S, 


j+1 


=  /XX. 


j+1 


(s)  [e  q  (x  . )  q  .  (s)  ] 


-1/4 


(1-P) r (s) ds 


Thus,  in  this  problem,  there  are  no  regular  integrals  to  consider. 

It  suffices  to  find  collocation  points  for  since  analogous 
results  for  Sj+1  will  follow  by  replacing  s  by  1-s  and  making  the 
appropriate  sign  changes  in  (4.37b).  Thus,  using  Laplace's  method 


we  approximate  (4.37a)  by 

.!  _-PjU-s)  2 


Sj  ~  /  e 
J  0 


~  -1/4 

[s‘q(Xj)qj (s) ]  (1-P) r j (s)ds  . 


For  taut  spline  approximations,  the  choice  of  given  by 

^  A  « 

(3.13b)  reduces  both  u^(s)  and  v^(s)  (cf.  (4.11))  to  0(h  /p^) 

and  it  is  still  reasonable  to  select  t.  according  to  (4.25),  i.e., 
•  * 

so  that  (4.24)  is  satisfied  with  and  given  by  (3.13b)  and 
(4.36b),  respectively.  We  continue  to  refer  to  this  choice  of 


tj  as  Method  1.  Likewise,  for  polynomial  approximations,  we  select 

tj  according  to  (4.31),  i.e.,  so  that  (4.30a)  is  satisfied,  and 

still  refer  to  this  as  Method  3.  Both  methods  reduce  S_.  (hence, 

S.  . )  to  at  least  0(ht.).  Since  t.  is  0(l/o.)  or  0(l/o.)  for 
J J  3  3  3 

Method  1  or  3,  respectively,  and  h/p^  and  h/p ^  are  both  0 ( /z) 
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(cf.  (3.13b)  and  (4.36b)),  we  would  expect  (cf.  (4.35a))  e(x^) 
to  be  at  most  0(h/£j  at  knots  away  from  boundary  layers.  The 


(5.1a) 


(5.1b) 


(5.2a) 


(5.2b) 


(5.  2c) 


computational  results  of  Section  5  support  this  conclusion. 

5.  Numerical  Results 

In  this  section,  we  apply  Methods  1,  2,  and  3  of  Section  4 

to  three  examples  having  known  exact  solutions.  Our  calculations 

are  performed  on  a  uniform  partition  with  spacing  h  and  the  adjoint 

★ 

tension  parameters  are  approximated  as 


Pj  =h 


lp(xk)A  -  q(xk)/p(xk)  |  ,  |up(Xj)/e  |  _>  |uq(Xj) 


/-yq(Xj )/2e 


|up(Xj)/e|  <  |uq(Xj) | 

j  =  1, 2, . . . ,N, 


where 


k  H 


j  -  1  if  up(Xj)/e  <  0 

j  if  yp (x j ) / z  >  0 


The  tension  parameters  and  collocation  points  t^  on  1^  are  chosen 
as  follows: 

Method  1:  Select  p^  according  to  (3.12)  and 

tj  =  min  {(1/Pj)  In  (l+p,/p*)  ,  (l-l//3)/2> 

Method  2:  Select  p.  according  to  (3.12)  and 

1  2^-12  P-i 

fcj  =  1  [1  '  ?"  COsh  (—  Sinh  T^1 

3  3 


Method  3:  Select  p.  =  0  and 
0  :  w(p*/2) 


t .  = 


J  PJ  1  +  /I-4u>('p*/2)/o* 

We  also  consider  "partial  tension"  methods  where  the  above 
rules  for  selecting  p^  and  tj  are  only  applied  on  subintervals 
containing  boundary  layers  and  collocation  is  performed  at  the 
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Gauss-Legendre  point  with  pj  =  0  on  all  other  subintervals.  These 

4 

methods  can  potentially  converge  as  0 (h  )  outside  of  boundary 


layer  subintervals  provided  that  the  numerical  solution  and  its 


derivative  have  been  computed  accurately  enough  at  the  ends  of 


subintervals  containing  boundary  layers.  Thus,  we  would  not  expect 


partial  tension  to  be  useful  with  Method  2  since  there  can  be  large 


errors  in  the  derivative  of  the  computed  solution  at  the  knots. 


We  denote  the  partial  tension  methods  that  use  either  Method  1 


or  3  within  boundary  layer  subintervals  as  either  Method  IP  or 


3P,  respectively.  In  order  to  automatically  locate  subintervals 


containing  boundary  layers  we  first  compute  a  preliminary  solution. 


c^,  j=0,l,...,N  using  either  Method  1  or  3  on  all  subintervals. 


Using  this  solution  we  calculate 


(5.3) 


=  |y  [f  (Xj)  -p(Xj)d^°-q(x..  )  c  ^ 0  ]  /2e  ,  j=l,2,...,N  . 


and  set  pj  =  0  and  tj  =  (l-l//3)/2  on  any  subinterval  having 


(5.4) 


_  »  «  _  *  » 


<  5  min  {l,yi  ,y2  , . . .  ,yN  >  , 


where  5  is  a  threshold  constant  which  we  normally  take  as  50.  The 


problem  is  then  re-solved  using  the  new  values  of  o  and  t..  This 

j  3 


procedure  is  somewhat  ad  hoc  and  has  not  been  totally  effective. 


e.  ecially  when  h  is  relatively  large  and  Cj0  and  dj^  are  inaccurate. 

_  I  I 

This  can  cause  errors  in  yj  which  can  lead  to  the  erroneous 


conclusion  that  a  subinterval  contains  a  boundary  layer  when  in  fact 


it  does  not,  and  vice  versa. 


Each  example  was  solved  for  various  values  of  e  and  h  =  (b-a)/N 


with  N  =  2  ,  k  =  2, 3,..., 7.  The  error  in  the  solution  and  its 


derivative  on  a  partition  with  spacing  h  are  measured  by 


(5.  5) 


xj£40 


yu)  (Xj)  -  y^11  (X.) 


i  =  0,1 


where  AQ  is  a  fixed  uniform  partition  that  is  specified  with  each 
example.  The  order  of  convergence  r^  is  computed  as 

(5.6)  ri  =  ln  (  I  I  e (i)  I  lh#  Aq/|  |e(i)  |  lh/2^  Aq>  /m  2  ,  i-0,1. 

All  calculations  were  performed  in  double  precision  on  an 
IBM  3033  computer.  Errors  that  are  less  than  5  x  10~^3  are  recorded 
as  zero  in  the  tables. 

Example  1; 

i  i  t  -v/2 

ey  +  ( ( 1+x) y)  =  e  x/  [  (1+x)  (3-x) +e/2] /2  ,  0  <  x  <  1  , 

(5.7) 

y (0)  =  0  ,  y(l)  =  e"1/2  -  e“7/3e 

The  exact  solution  of  this  problem  is 

2 

y(x)  =  exp  (-x/2)  -  exp  [-x(x  +3x+3)/3e]  . 

There  is  a  boundary  layer  of  width  0(e)  near  x  =  0. 

In  order  to  demonstrate  how  poorly  collocation  at  the  Gauss- 

Legendre  points  with  cubic  polynomials  can  behave  on  singularly- 

-4 

perturbed  problems,  we  solved  (5.7)  with  e  =  10  and  h  =  1/8  by 

this  method  and  plotted  the  computed  solution  in  Figure  2.  It 

bears  little  resemblance  to  the  exact  solution  which  is  essentially 
-x/2  -3 

e  for  x  >  10  .  Pointwise  the  numerical  solution  approximately 

lies  on  the  straight  line  joining  the  two  boundary  values  y(0) 
and  y (1)  . 

We  solved  (5.7)  for  e  =  10  1,  i  =  2, 4, 6, 8  using  Methods  1,  2, 

I 

3,  and  IP.  The  errors  | je] |h  ^  and  I le  11^  a  on  the  Partition 

Ag  =  (1/8,  1/4,  3/8,...,  7/8}  are  presented  in  Tables  1  and  2, 

—2  -4  —8 

respectively,  for  z  =  10  ,10  ,  and  10  .  (The  results  for 

e  =  10  were  essentially  the  same  as  those  for  10  . )  | |e| ], 

h'u0 

is  also  plotted  as  a  function  of  1/h  in  Figure  3  for  Methods  1,  2, 
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and  3  and  e  =  lO”1,  i=2,4,6,8.  Tables  1  and  2  indicate  tnat  when 

2 

e/h  <<  1,  Methods  1,  2,  and  3  are  converging  as  0(h),  0(h  ),  and 

3  *  2 

0(h),  respectively,  and  that  ||e  ||.  .  is  converging  as  0 (h  /e) 

n,  Ao 

_2 

for  Method  2.  For  larger  values  of  e,  e.g.  e  =  10  ,  the  order 
of  convergence  can  be  seen  to  increase  as  h  decreases  (e/h  increases) 
and  the  collocation  points  move  closer  to  the  Gauss-Legendre  points. 
Partial  tension  with  Method  IP  yields  a  dramatic  improvement  in 
the  results  obtained  by  Method  1. 

In  order  to  provide  some  indication  of  how  Methods  1,  2,  and 
3  behave  on  subintervals  containing  boundary  layers  and  between 
knot  points,  we  have  plotted  their  computed  solutions  y^(x)  on 

I 

0  <  x  <  2h  (Figure  4)  and  their  errors  e(x)  and  e  (x)  on  0  <_  x  <_  1 

-4 

(Figures  5,6,7)  for  e  =  10  and  h  =  1/8.  The  error  in  the  Method  2 
solution  shown  m  Figure  4  is  less  than  3.2  x  10  for  all  x  e  I^u^- 
Method  1  yields  poor  accuracy  for  xel.^  but  it  does  predict  the 
correct  width  of  the  boundary  layer.  Method  3  dissipates  the  boundary 
layer;  however,  the  dissipation  is  largely  confined  the  subinterval 
1^  containing  the  boundary  layer.  Figures  5-7  show  that  all  three 

I 

methods  have  spurious  internal  boundary  layer  jumps  in  y^(x)  at 
the  ends  of  each  subinterval  and  that  Method  2,  because  of  the 

I 

singular  behavior  of  yh  in  e  at  the  knots,  has  spurious  jumps  in 
yh(x)  itself.  The  jumps  in  y^(x)  are  small  and  one  is  not  normally 
interested  in  the  solution  at  points  other  than  the  knot  points; 
however,  the  results  indicate  that  some  caution  is  needed  when 
using  (2.7-12)  to  interpolate  the  solution  between  knot  points. 

Figures  5-7  further  indicate  that  the  pointwise  error  is 
largest  at  x^  =  h  and  decreases  at  knots  that  are  farther  frcm 
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(5.8) 


the  boundary  layer.  This  is  generally  true  for  other  values  of  h 

as  well  and,  thus,  we  have  tabulated  |e(h) |  for  Methods  1,  2,  and 

3  with  £  =  10  1,  i  =  2, 4, 6, 8  and  N  =  21,  i  =  2,3,.. .,7,  in  Table  3. 

The  results  for  Methods  2  and  3  indicate  that  |e(h)  j  cannot  be 

2  3 

reduced  below  0(e)  until  e/h  and  e/h  ,  respectively,  are  sufficiently 
small. 

Example  2; 

2  '  '  2  2  2  2 
e  y  -  (x  +e)y  =  -  (x  +e)  (1+sin  irx)  -  e  ir  sin  irx,  0  <  x  <  1  , 

y(0)  =  y (1)  =  0  . 

The  exact  solution  of  this  example  is 


y(x) 


1  +  sin  irx 


- - -  {  (l-e“1/2e)W(x//e) e"x  /2e 

erf  (/i/e) 


+  [l-e"1/2eW(/T7e) ]e"(1_x2)/2e}  , 

where 

z2 

W  ( z)  =  e  erfc(z)  . 

»  *  2 

(Note  that  y  is  multiplied  by  e  instead  of  e  for  notational 

simplicity.)  The  exact  solution  of  (5.8)  has  a  boundary  layer 

of  width  0(/e)  near  x  =  0  and  one  of  width  0(e)  near  x  =  1.  This 

problem  does  not  satisfy  the  assumptions  of  Problem  2  since 
2 

q(x)  =  -(x  +e)  cannot  be  bounded  away  from  zero  at  x  =  0  as  e  •*  0; 
thus,  x  =  0  is  a  turning  point. 

We  computed  solutions  by  Methods  1  and  3  for  e  =  10-1,  i=2,4,6,8. 
Solutions  were  also  calculated  using  the  collocation  points  of 
Method  3  (cf.  (5.2c))  but  with  splines  under  tension  instead  of 
polynomial  splines,  and  these  are  denoted  as  Method  3T.  The  errors 
|  |e|  a  and  !!e  |  |j^  ^  are  presented  in  Tables  4  and  5,  respectively, 

for  the  partition  AQ  =  (2/8,  1/4,  3/8  , .  . .  ,  7/8  }  .  Partial  Tension 


28 


1 


solutions  using  Methods  1  and  3T  were  also  calculated  and  the 

results  were  marginally  more  accurate  than  those  of  Method  1. 

Tables  4  and  5  indicate  that  when  e/h  <<  1,  |  |e|  |,  .  and 

1 1  1  h,AQ 

i  2 

||e  | |.  A  are  0{eh)  and  0(h),  respectively,  for  Method  1  and 
n, 

3  2 

0(eh  )  and  0(h  ),  respectively,  for  Method  3T.  The  small  errors 
in  Method  3  make  it  difficult  to  estimate  the  order  of  convergence. 

The  largest  error  on  the  partition  A  =  {0,h, 2h, . . . ,Nh  =  1} 
used  for  the  computation  was  once  again  at  the  end  of  a  subinterval 
containing  a  boundary  layer,  i.e.,  either  at  =  h  or  xN-1  =  (N-l)h. 
To  indicate  how  this  error  behaves,  we  present  results  for  | |e| |.  , 

in  Table  6.  For  small  values  of  z/h  we  see  that  the  polynomial 
solution  (Method  3)  is  not  converging  in  h  and  that  this  situation 
is  remedied  by  using  taut  splines  (Method  3T) .  Although  we  have 
not  done  so,  we  suspect  that  it  would  have  been  sufficient  to 
use  taut  splines  only  within  subintervals  containing  boundary 
layers.  Furthermore,  Methods  1  and  3  do  not  appear  to  be  uniformly 
convergent  in  h  for  all  e.  All  methods  are  converging  as  0(e) 
which  accounts  for  the  remarkable  accuracy  when  e  is  small. 

Example  3:  (cf.  Hemker  [15]) 

i  i  i  2 

ey  +  xy  =  -  eir  cos  ttx  -  ttx  sin  ttx  ,  -1  <  x  <  1  , 

(5.9) 


y (-1)  =  -2  ,  y(l)  =  0  . 


The  exact  solution  of  this  example  is 

y(x)  =  cos  irx  +  erf  (x//Te)  /erf  ( l//2z )  . 

The  problem  has  a  turning  point  at  x  =  0  and  the  exact  solution 
features  an  interior  or  "shock"  layer  there. 

Solutions  were  calculated  by  Methods  1,  2,  3,  and  IP  for 
£  “  10  1,  i=2,4,6,8.  The  errors  |  |e|  |,  .  and  |  j  e  |  | .  .  on  the 

'  h, Aq  1 1  '  h, Aq 
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partition  AQ  =  (-3/4,  -1/2,  -1/4,  1/4,  1/2,  3/4}  are  presented 

-2  -4  -8 

in  Tables  7  and  8,  respectively  for  e  =  10  ,10  ,  and  10 

(The  results  for  e=  10-6  were  essentially  the  same  as  those  for 
_8 

e  =  10  .)  The  results  for  ||e||.  A  are  indicating  the  same 

orders  of  convergence  as  found  in  Example  1;  however,  for  Methods  1 

and  3,  [  |e  |L  *  is  much  more  accurate  than  the  corresponding 

n,aQ 


values  for  example  1. 

In  order  to  include  the  behavior  of  the  solution  in  the  turn¬ 


ing  point  region,  we  have  tabulated 


1  h,  A 


on  the  partition 


A  =  {-1, -1+h, . . . , -1+Nh=l}  used  for  the  calculation  in  Table  9. 


Note  that 


h.  A, 


h,  A 


for  Method  1,  so  the  maximum  error 


is  not  in  the  turning  point  region.  Methods  2  and  3  are  both 
exhibiting  regions  of  non-uniform  convergence  in  h. 

6.  Discussion  and  Conclusions 

Based  on  the  results  of  Sections  4  and  5,  we  conclude  that 
splines  under  tension  are  most  suitable  in  regions  containing 
boundary  and/or  interior  layers  and  that  collocation  with  piecewise 
polynomials  are  superior  elsewhere.  In  particular,  when  e/h  <<  1 
Method  3  provides  an  approximation  that  converges  outside  of  boundary 
layer  subintervals  as  O(h^)  when  p(x)  ^  0  and  at  least  0(/£h) 
when  p(x)  ^  0  on  [a,b] .  Hemker  [15]  and  de  Groen  and  Hemker  [10] 
reached  similar  conclusions  with  their  exponentially  fitted  Galerkin 
methods. 

4 

Partial  tension  can  converge  as  0 ( h  )  outside  of  boundary 
layer  subintervals,  but  either  requires  a  knowledge  of  boundary 
layer  locations  or  a  preliminary  solution  to  automatically  locate 
them.  The  latter  procedure  may  be  useful  for  nonlinear  problems 
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where  it  is  necessary  to  solve  a  sequence  of  linear  problems  to 
find  the  solution;  however,  for  linear  problems  it  does  not  seem 
to  warrant  the  extra  computational  effort  merely  to  obtain  one 
order  of  accuracy  more  than  that  available  by  Method  3. 

The  use  of  "one  sided  splines  under  tension,"  i.e.,  the 
selection  of  basis  functions  that  satisfy 

I  III 

(6.1)  (n  +  pn)  =o  ,  o  <  s  <  l  , 

instead  of  (2.14)  would  undoubtably  improve  the  results  on  sub¬ 
intervals  where  p(x)  ^0.  A  basis  for  these  approximations  would 
contain  either  the  exponential  exp(-ps)  when  p  >  0  or  exp(- | p | (1-s) ) 
when  p  <  0,  and  not  both  as  in  the  current  case.  This  would  more 
accurately  represent  the  exact  solution  of  problems  where  p(x)  ^0 
on  [a,b].  The  results  could  possibly  be  further  improved  by  not 
restricting  the  collocation  points  to  be  placed  symmetrically 
on  each  subinterval  and  by  using  non-uniform  partitions.  Each 
of  these  potential  improvements  is  currently  under  investigation. 
Extensions  of  our  methods  to  higher  order  scalar  and  vector  systems 
of  two-point  boundary  value  problems  as  well  as  second  order 
parabolic  partial  differential  equations  are  also  being  studied. 
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Table  1:  Error  and  Order  of  Converaence  for  Example  1  Measured  on 
-  U/8,  1/4,  3/8,...,  7/8} . 
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Table  2:  Error  and  Order  of  Convergence  in  the  Derivative  of  the 

Solution  of  Example  1  Measured  on  AQ  =  (1/8,  1/4,  3/8,..., 7/8} 
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Table  3:  Error  at  the  Knot  Point  x^  for  Example  1. 
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Table  4:  Error  and  Order  of  Convergence  for  Example  2  Measured  on 
AQ  =  {1/3,  1/4,  3/3,...,  7/8}. 
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Table  7:  Error  and  Order  of  Convergence  for  Example  3  Measured  on 
AQ  =  {-3/4,  -1/2,  -1/4,  1/4,  1/2,  3/4}. 
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Table  8:  Error  and  the  Order  of  Convergence  in  the  Derivative  of  the 
Solution  of  Example  3  Measured  on 
AQ  =  {-3/4,  -1/2,  -1/4,  1/4,  1/2,  3/4}. 
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Table  9:  Error  and  Order  of  Convergence  for  Example  3  Measured  on 
A  «  {-1,  -1+h,  -l+2h, . . . ,  -1+Nh  =  1}. 


Figure  Captions 


Figure  1:  Canonical  basis  functions  nQ(s,p)  and  n1(sfp)  on  -1  <  s  <  1 
for  p  =  0.01  (Figure  la)  and  p  =  10  (Figure  lb). 

Figure  2:  Solution  of  Example  1  using  cubic  polynomials  and 
collocation  at  the  Gauss-Legendre  points. 

Figure  3:  Error  |je(h)||.  for  Example  1  using  Methods  1,  2,  and 

3  and  measured  on  the  partition  AQ  =  (1/8,  1/4,  3/8,...,  7/8} 
Figure  4:  Solutions  of  Example  1  using  Methods  1,  2,  and  3  on 
0  £  x  £  1/4  for  e  =  10  ^  and  h  =  1/8. 

I 

Figure  5:  Error  e(x)  and  its  derivative  e  (x)  in  the  solution  of 
Example  1  by  Method  1  for  e  =  10  4  and  h  =  1/8.  (Note: 
e' (0)  =  -5.76  x  103. ) 

I 

Figure  6:  Error  e(x)  and  its  derivative  e  (x)  in  the  solution  of 
Example  1  by  Method  2  for  e  =  10  4  and  h  *  1/8. 

I 

Figure  7:  Error  e(x)  and  its  derivative  e  (x)  in  the  solution  of 
Example  1  by  Method  3  for  e  =  10  4  and  h  =  1/8.  (Note: 
e' (0)  =  1.0  x  104.) 
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